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ABSTRACT 

The main issue of the pulsar magnetosphere is how the rotation power is converted into 
both particle beams which causes pulsed emissions, and a highly relativistic wind of 
electron-positron plasmas which forms surro unding nebulae shining in X-rays and TeV 



gamma-rays. As a sequel of the first paper (jWada fc Shibatall2007() . we carried out a 
three dimensional particle simulation for the axisymmetric global magnetosphere. We 
present the results of additional calculations, which are higher resolution model and 
higher pair creation rate cases, and a detailed analysis for the solution. We confined to 
demonstrate the cases of low pair creation rate, i.e., the magnetic field is fixed dipole. 
The radiation drag of the plasma is taken in a form with the curvature radius along the 
dipole magnetic field. The electrostatic in teractions are calcula ted by a programmable 



special purpose computer, GRAPE-DR (jMakino et al.ll2007l ). Once pair creation is 



onset in the outer gaps, the both signed particles begin to drift across the closed 
magnetic field due to radiation drag, and they create outflow. Eventually, the steady 
magnetosphere has outer gaps, both signed outflow of plasma and a region in which 
the electric field is dominant extending from the equator. In the steady state, the 
magnetic field made by magnetospheric current is comparable to the dipole magnetic 
field outside of several light radii from the star. In much more pair creation rate model, 
the effect of modification of the magnetic field will bring about modification of the 
outflow of the plasma, requiring further study with higher pair creation rate model in 
a subsequent paper. 

Key words: pulsars: general - magnetic fields - plasmas. 



1 INTRODUCTION 



Rota tion powered pulsars are bright sources in the X-ray sky (e.g., Becker fc Triiempe J 19971 : iPossenti et al iKargaltsev fc Pavlov! 



2005 ) and the gamma-ray sky re.g.. lNel et allll996l : de Jager fc Diannati-AtailboOSl l as well as in the radio (e.g. 



Hobbs et al 



200J), and they are identified as magnetized rotating neutron stars. The pulsed emissions originate in the accelerated particles 



in the magnetosphere. In addition, as a persistent source of hi gh energy photons, we observe neb ulae around young pulsars 
with synchrotron X-ray radiation ( Kargaltsev fc Pavl ov 20081'). inverse Compt on TeV radiation ( de Jager fc Diannati-Atai 
20081 ). and thermal radiations from the neutron star ( Becker fc Triiemper 1997 ). The emission from the nebula is caused by 



termination shock of t he pulsar wind, which is believed to be relativis tic outflow of magnetized pair plasmas carrying out most 
of the rotation power (|Rees fc Gunnlll974l : iKennel fc Coronitil Il984al lb[) . However, the nature of the pulsa r wmd and the radi- 
ation beam from the pulsar magnetosphere have remained a longstanding problem. As a primitive model. iGoldreich fc Julian 
19691 . hereafter GJ, suggested an axisymmetric steady model of pulsar magnetosphere in which the neutron star is assumed 



to be a perfect conductor surrounded everywhere by charge-separated force-free plasmas. The model anticipates that the 
plasma holds corotation with the star in the light cylinder although this corotational magnetosphere can be violated beyond 
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the light cyhnder. Then the GJ model assumes that the plasma is replaced by the acceleration along the opened magnetic 
field line beyond the light cylinder. On the basis of the copious plasma in the magnetosp here, the global re lat ivistic mag- 
netodynamics models are ap plied to help us understand acceleration of the pulsar wind ( McKinnev 2006al lbl: iKomissarov 



20061 : 



Bucciantini et alj 



20061) and possible dissipation such as in magnetic neutral sheets may contribute to high energy 



pulsed radiation (jXirk et'liL 2002 '). Against these studies, several studies showed that an aligned rotating neutron star should 
consist of the electrostatic charge clouds surrounded by the vac uum gaps without pul sar wind and therefore the steady 
solutions with outflow of t he pla sma are not self-consistent (See. iMichel fc Smithll200ll ). About the electrostatic solution, 
Krause-Polstorff fc Michel (|l985al lbl. hereafter KM) performed a particle simulation for the axisymmetric model and showed 



formation of the gaps around the null surface with a rotating equatorial disc and polar domes. This model suggests a quiet 
electrosphere when the plasma is extracted only from the ste llar surface, but the vacuum gap in the middle latitudes looked 
to be unstable against the pair creation cascade. Successively, ISmith. Michel, fc Thackeii (|2001f ) reconstructed the same result 
with a higher accuracy and even if the pair plasma were generated in the gaps at the middle latitudes, the quiet solution 
without pulsar wind was identified. In their models, the quiet solution is obtained that the particles are put at equilibrium 
positions along the dipole magnetic field line and this picture would be valid if the plasma is set down fully inside the light 
cylinder and the magnetic field is larger than the electric field. In other words, the effects of the inertia of the c o-rotat ional 
plasma around the light cylinder is not considered. The same result is obtained by Petri. Hevvaerts. fc Bonazzolal ( 2002bh who 
concluded a given total charge of the system determined by the net outgoing charged flux and the potential configuration is 
unfavorable for the particles to escape to infinity such as the active solution obtained by MHD approximation. However, the 
quiet solution should be modified to concern the leaking of the plasma at which the electro static plasma is in the vicinity of 
the light cylinder. As a possibility of radial outflow from the equatorial disc on the equator. ISpitkovskvl (|2004 ) discovered the 
leakin g mechanism of the equatorial disc caused by the diocotron instability with numerical simulation, which is also pointed 
out bv lPetri et al.l (|2002al ) with linear theory. The leaking of the disc should be significant to make the current in the vacuum 
region outside of the quiet electrosphere. But the time scale of the diffusion of the plasma in the vicinity of the light cylinder is 
much larger than the rotation period and therefore further study is needed to follow the global current structure with a much 
longer time scale such as a global numerical simulation with much more compu tational time with relat ivistic kinematics. 

Contrary to the previous quiet solutions, we reported an active solution ( Wada fc Shibatal 2007 . hereafter, WS) with 
coexisting pulsar wind and outer gaps in which the effect of the inertia and radiation drag for the plasma is included by 
3-dimensional particle simulation with special relativistic regime, and the simulation starts from the quiet solution as KM. 
Once pair plasmas are generated in the vacuum gaps around the quiet electrosphere, the positive charged outflow is formed 
near the equator and the negative charged outflow is formed from the polar region of the star. The outflow of plasma 
maintains the charge deficiency in the outer gaps. For the balancing of loss and the supply of particles from the outer gaps, 
the magnetosphere eventually settles in a steady state with pulsar wind. But there remains some problems to solve. At first, 
although our simulation showed that the wind coexists with the pair creating regions in the pulsar outer magnetosphere, the 
solution is only obtained by a lowest case of pair creation rate, where we carried out the simulation in which the current 
density around the star is less than the one of GJ model, and therefore magnetic field by the magneto-spheric current is 
negligibly-small in the light cylinder. If the pair creation rate taken to be large gradually in our simulation, the gaps tend to 
reduce and decrease the generation rate of the plasma and therefore the active structure might be disappeared. Although the 
drift motion of plasma caused by radiation drag force works as a leaking mechanism in closed magnetic field and is favorable 
to maintain the outer gaps, recalculation should be needed whether our active solution is maintained in higher pair creation 
rate models. Secondly, the obtained structure of the outer gaps contains some artifact due to the accuracy. The detailed flow 
pattern of plasma such as the poloidal current loop inside of the light cylinder was not clear because of low resolution of our 
previous work. We will simulate the lowest pair creation rate model again with higher resolution in which the unit charge 
of the particle is taken to be 1/10. At last, our particle simulation has numerical error, which is defined as the resolution of 
charge of simulation particle and the injection of the plasma from the inner boundary. In our previous work (WS), the polar 
potential drop is maintained in three times larger fraction than artifact level, present higher resolution model will decrease 
the error and could identify the existence. However the particle has large inertia length which is about 3% of the stellar radius 
and therefore we denote that our simulation could not resolve the size of polar cap accelerator. 

Meanwhile the development of global magnetosphere, some phenomen ological models which concerns the local structure 
for the magnetospheric activit i es have been developed: the p olar cap models ( Ruderman fc Sutherland! 19751: 



19791: IZhang fc Hardinjboool: [Harding fc Muslimov 



: tne p olar cap moaels l|Jri,uaerman &: Sutnerianclllb)Ya:lArons fc Scharlemann 
2002!) ■ the outer gap models (lHollowavlll973l: Cheng. Ho. fc RudermanI 



1986al lbl: iRomanil 1 19961; Zhang fc Chenj 11997*: 'H irotani fc Shibatal [l999l : iTakata. Shibata. Hirotani. fc Chanel booel 'l and the 
slot gap model ( Muslimov fc Hardinj 200S . i2004il have been developed to explain the observed pulsed light curve and spec- 
trum. The critical issue for these models is how the electric field along the magnetic field, By, is maintained locally although 
the existence of the gap is implicit in these models. Although the adequacy of these gap models is achieved statistically by 
many further observations, but these models usually contain many degrees of freedom, for example, which are the surface 
magnetic field intensity of the pulsar, the inclination angle between magnetic axis and rotational axis, the viewing angle from 
the gaps and so on. Recently, the radiation from the pulsar magnetosphere with high-energy 7-ray bands has been observed 
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by the Fermi 7-ray telescope, and the 38 gamma-ray pulsars have been discovered ( Abdo et alj l2009al lbl. boiol ). including 
21 radio-loud and 17 radio-quiet. These observations could discriminate between the emission models located on the polar 
region or outer region, i.g., polar cap, slot gap, and outer gap. The Fermi telescope has also measured the spectral properties 
above 10 GeV with a better sensitivity than EGRET. It was found that the spectral shape of 7-ray emissions from the Vela 
pulsar is well fitted with a power low (photon index T ~ 1.5) plus exponential cut-off (- Ecut ~ 3 GeV) mod el. The discovered 
exponential cut-off feature predicts that the emissions from the outer magnetosphere (jAbdo et al.ll2009d ) are more favored 
than the polar cap region jPaughertv fc Harding||l996l ). which predicts a super exponential cut-off with the magnetic pair- 
creation. Furthermore, the detection of the radiation a bove 25 GeV ban ds associated with the Crab pulsar has also predicted 
the high-energy emission in the outer magnetosphere ( Aliu et al. 20081 ) with MAGIC (Major Atmospheric Gamma Imaging 
Cherenkov) telescope. 

In this paper, we develop our global model such as linking the outer gap with the pulsar wind. In order to construct a 
solution, we utilize particle simulation in which particle motion and electric field are alternatively solved in the same way 
as the particle-in-cell (PIC) methods but with the magnetic field given with dipole in our present model. We simulate the 
generation of the plasma by pair creation in the gap and then obtain a steady state of the axisymmetric magnetosphere by 
solving the equation of motion including drag force due to curvature radiation and the electromagnetic fields. Since the motion 
of the individual particles is tracked, we can simulate any kind of drift motions for plasma, arising from the radiation drag 
force, the centrifugal force and the gradient of the magnetic field; f^^^ x B drift, E x B drift and magnetic gradient drift. 
Here, we provide results of additional calculations and the detailed analysis for our steady solutions. We shall show that the 
trans-field drift by radiation drag plays an important role and eventually the electric field dominant region expanding to both 
sides of the equatorial plane makes a global poloidal current loop. But this paper considers cases of low rate of pair creation, 
and therefore modification of the original dipole field would be ignored as well as in our previous work. At first, the previous 
result (WS) of our steady solution was reconstructed with 10 times number of the particles and models with several pair 
creation rates were carried out. Our methods of calculations for the particle simulation are given in section (2] and the result 
is presented in section [3] and section [4] is for discussion. 



2 NUMERICAL METHOD 
2.1 Outline 

We solve the motion of the particles and the electromagnetic field alternatively so as to obtain a self-consistent steady 
structure of the electromagnetic field and particle distribution. This method is similar to the well-established particle-in-cell 
(PIC) method. However, to obtain a steady solution, we use static solutions of the Maxwell equations, i.e., we omit the effects 
of time variation of the field ( dE/dt — dB/ dt = 0). This enabl e us to use the programmable special purpose computer for 



astronomical N-body problem; GRAPE-DR (jMakino et al.l 120071 ). which calculates Coulomb interaction very rapidly with a 



high degree of precision at the position of the plasma. 

In our simulation, generation of the particle is considered, which is not common in PIC methods. There are two cases: 
(1) electrons or protons on the stellar surface pop into the magnetosphere due to the unipolar induced electric field of the 
star, (2) electron-positron pairs are created due to photon-photon collision or magnetic pair creation in the magnetosphere. 
For the simplicity, the mass of both signed simulation particles is taken to be the same value, that is, the proton and positron 
are not distinguished, and the type of pair creation is not distinguished, that is, the detailed process of pair creation is not 
considered. We assume the pulsar to be a rotating spherical conductor in which magnetic field is uniform with the magnetized 
axis being parallel to the rotation axis. In this paper, because we consider the case of low rate of pair creation, which has 
much less current density than GJ model inside light cylinder, the modification of the dipole magnetic field near the star by 
the magneto-spheric current would be trivial. For this reason, the magnetic field outside of the star is assumed to be dipole, 
i.e. 

B = ^(2 cos ee,. singes), (1) 

where 6 is the colatitude, r is the distance from the center of the star, jj, — BqR^ /2 is the dipole magnetic moment, Bq 
is magnetic field intensity at the poles, R is the stellar radius and and ee are the unit vector along r and 9 directions, 
respectively. These assumption gives the electric potential on the stellar surface, 

, , sin" 9 

= ^ h constant, (2) 

where SI is the angular velocity of the star. This provides the inner boundary condition for the Poisson equation for the electric 
potential. 

If the outside of the star is a vacuum, the electrostatic potential given by the solution of Laplace equation V^(^ = is 

(pv = 4>i + <t>i, (3) 
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where 

(4) 



3cr 



r 

2 

(Scos^f-l), (5) 



and Qaya = a2iin/(3c), where a is the non-dimensional net charge of the system. If we foUow the Jackson's Gedanken 
experiment ( Jacksonlll976l ) , we choi 
the vacuum electric field becomes 



periment ( Jacksonllf 976l ). we choose a = 1 as an initial condition, which corresponds to the +10 model for KM. Consequently, 

E^r) = -V0V = Ei{r) + Ei{r), (6) 



where 



EA{r)^^i^-—{'icos'e~l)er-- r ee, (7) 



E,{r)^^er. (8) 



r 



This vacuum solution has the surface charge density where 



^ 2 Q\ I Qsys 



Given Bo, R and the light speed c, all quantities in our simulation are normalized. In the following, the barred quantities 
indicate that they are normalized values. For example, electric charge is measured in units of BqR? , and Qsys = Qsys/ B^R? . 
In the following subsections, we describe the detailed method of our simulation. 



2.2 Method of Solution for the Electric field 

The space charge density in the magnetosphere is represented by limited number of simulation particles, i.e.. 

n 



(10) 



In simulation, we are concerned with the super particles which have the same mass and opposite signed charge with the same 
absolute value of the charge and therefore the difference in the mass of ion and that of positron is not considered. As we 
introduce later, the normalized value of the mass and charge for the super particle is represented by m and q, respectively. 
The electric potential in the magnetosphere is determined by Poisson's equation, i.e., the solution for electric potential is 
given by the superposition of the vacuum component ([3]) and space charge component, which is calculate by — V'^(/)m ~ Pm, 
with the boundary condition <j)n-i{R) — 0, and we have 

1 



for the electric potential. 



for the electric field, and 



R/r 



{R/n 



r, R r-iR/nfn 



1-^U 



{R/n 



R^ 



1-^ 

ri 



(11) 



(12) 



(13) 



ji-KR\R-ri\'i 47ri?2 

for the surface charge density. Thus, the solution of the electric field with the boundary condition ^ is given in the forms, 

= (?!)m + ^v, (14) 
E = En, + Bv, (15) 
O- = CTm + CTv. (16) 



2.3 Popping of particles from the star surface 

GJ pointed out that the rotation-induced electric field pulls charged particles from the stellar surface against the gravitational 
force. In the following, we describe how this process is realized in our simulation. 

In vacuum conditions, the scalar product of the electric field and the magnetic field on the stellar surface is given by 

i5.B = ^^cose(|-cos^e). (17) 

The sign of the electric field along the magnetic field (-Ey) changes on the surfaces cos^ = ±\/a/3, i.e., the negative charged 



A Particle Simulation for the Pulsar Magnetosphere: 11 5 



particles pop out from the polar region and positive charged particles pop out from the equatorial region. Once the charge is 
emitted from the stellar s urface, they e xperience _E|| which changes sign on r = i? + \/2i/a cos 6 for electron or on 61 = 7r/2 
for positron and ion (See, Jackson 19761 ) . The electrons are emitted from the lower colatitude region are accelerated outward 
along the magnetic field line at first. For electron, after passing over the surface at which the minimum of the potential energy 
is given, the field aligned electric field decelerate the particles and they are refiected toward the stellar surface again. Then 
they accumulate near the surface and makes the polar dome of electrons because of the radiation drag force by the curvature 
radiation for the particle. At the same time, positrons or protons are emitted from the higher colatitudes accumulating on the 
latter surface, i.e., the equatorial plane, and making the equatorial positively-charged torus. Thus, the vacuum electric field 
tends to be screened out by the particles in the magnetosphere. 

If the surface electric field along the magnetic field is shielded, the surface charge density becomes, 

3 Bq^R . 2. 

o-G,i = T, sm e. (18) 

87r c 

This surface charge density is resulted in the kink of the magnetic field from the uniform one inside to dipole one outside 
under the ideal MHD condition. Although the surface charge density appearing on the stellar surface is cr, only the excess 
charge a — ctgj is emitted from the stellar surface and therefore we expect that the surface charge density become ctgj in 
the steady condition. Thus, we replace a — ctgj by the movable simulation particle. For simplicity, the work function of the 
particles on the stellar surface is not considered, that is the free emission of particles is assumed. 



2.4 Parameter setting 

In our simulation, we use the super particle with artificially-enlarged mass and charge of g = 10~^,m = 10"^° or g = 
10~®,m = 10~^^. The number of particles representing the charge pG3,po\eR? is about 10* particles for q = 10~^ model, which 
is same setting with WS, and 10^ particles for q = 10~^ model, respectively. Actual numbers of the particle is increased by 
pair creation, the total number of particles in the simulation domain jumps tenfold for them and therefore it needs heavy 
number of calculation related to square of the number of particle for Coulomb interaction in one step. Thus, we gain the 
privilege of using the GRAPE-DR, which is a programmable special purpose computer for Astronomical N-body problems, 
to calculate interaction between the particles. 

In a realistic pulsar magnetosphere, the minimum scale of the motion of the plasma is gyro radius, which is rg = 
1.7 X 10^^r377cm for an electron or a positron, where 7 is Lorentz factor of the particle rg = {r/Rf, and 77 — 7/10^ 
respectively. We use super particles with a larger mass-charge ratio than the real one. We take H = 0.2 and fh/q = 10~^. 
For the definition of the numerical parameters, a typical gyro radius of the simulation particle on the pole is set to be 10~^ 
stellar radius, and therefore the corresponding time step is set in At = 10^^. Since the valid time step of particle for the gyro 
motion increases proportionately to the distance cubed from the star ~ f'^, we use the individual time step for each simulation 
particle adapted to the position to reduce the total numbers of the integration time. As we introduce the size of numerical 
domain later, we have to follow these particles over several dynamical time length of the simulation domain until a steady 
condition is achieved. It takes about 6 rotation periods of the star, i.e., tsim ~ 200. 

For the unipolar induced electric potential by the rotational magnetized conductor, the available maximum energy of the 
particle is estimated by the open field line voltage 

0.«^^^1.6xlO-(^)"\olt, (19) 

which gives maximum Lorentz factor of 

7_ = ^ = 3.2 X 10« f-^) (20) 
mc* \ 0.2 sec/ 

for electron and positron. Note that if the Lorentz factor of a particle becomes 7max, the gyro radius in the vicinity of light 
cylinder becomes the same order of the light radius. In other words, as expected, the localized acceleration region has much 
smaller potential drop than (f>as, the gyro radius of the particles is much smaller than the light radius. In our simulation, the 
gyro radius fg of all accelerated super particles is guaranteed smaller than the size of the light radius; typically fg — 10^'^-Ri 
for an accelerated particle corresponds to the Lorentz factor 0.l7max at the light radius. 



2.5 Inner and Outer boundaries 

For the injection of the particles from the stellar surface, the surface charge (a — ctgj) AS is replaced by the simulation particles, 
where AS is the surface element of the stellar surface. If the particles are set on the stellar surface, they are pulled back by the 
mirror charge, i.e., the second term of (|12[1. If this attractive field larger than the induced electric field, it prevents extraction 
of particles from the stellar surface as an artificial work function. Then the assumption of free emission of the particle from 
the stellar surface is not achieved. We put the particles above altitude of he from the stellar surface to avoid this artificial 



6 T. Wada and S. Shibata 



attractive force. We chose he — 0.04 for q = 10^^ model, he = 0.013 for q — 10^^ model, respectively. If the particles go back 
into the sphere with radius r — R + hc, the particles are deleted in the numerical region. It is notable that the numerical thin 
vacuum layer does not change the result at all. 

To save the total number of simulation particles, we remove the particles beyond a sphere with the radius -Rob and the 
removed charge is reduced from Qsys. The size of outer boundary is taken to be larger than the one of the accelerating region 
(~ Ri), that is, we chose _Rob = 10-Ri. In our simulation, the removed particles beyond the outer boundary have positive total 
energy (kinetic plus potential), and we consider they would be outflow. To examine the effect of the size of the outer boundary, 
the larger outer boundary models, which are -Rob = 20-Ri, 30-Ri are investigated, then we confirmed that the structure of the 
magnetosphere inside the sphere with radius 8Ri is not affected by the sizes of the outer boundaries. Thus, we use -Rob = 10-Ri 
as our standard model. 



2.6 Particle motion 

The equation of motion of the particles is given by 

rn.c^ = [S(rO + X B(rO] + /..d,,, (21) 

where rrii, qi, Vi and /3j are the mass, the charge, the position and the velocity of the i-th particle and f-^^^i i indicates the 
drag force due to curvature radiation of the accelerated particles. In this paper, the magnetic field is not deformed by the 
magneto-spheric current. The radiation drag force of the particles is given 

where 7?c is the curvature radius along a dipole magnetic field line at the position of the particle. Substituting vj'^ jr''' — 
l/''cq = const, the curvature radius is given by 

^ '^SihT? ^ = 2-r/re, ' 

where r^q is axial distance of the point at the intersection of dipole magnetic field line and the equatorial plane. Because 
1 < / < 4 with < r/rcq < 1, we take Rc = 4r/(3sin9) approximately. 

If the Lorentz factor of the particle increases, (|22|) become comparable to the Lorentz force, and such particle has f^^^ x B 
drift motion crossing the dipole magnetic surface. The critical Lorentz factor for this effect is given by 

For a realistic pulsar magnetosphere, it corresponds to 7d — 3.6 x IO^'Pq 2^''''(-Rc/-Ri)^'''^. Considering the dependency of 
the distance for -Rc, -Rc ~ -Ri in the vicinity of the light cylinder. For electron, 7d/7max = 0.11 with P — 0.2 sec and 
fi — 10^" gauss cm'', and thus the particles accelerated with 11% of the open field line voltage suffer from the such drift motion. 
On the other hand, for our simulation particles, 7d,sim/7max,sim ~ 0.0055 for q = 10~^ model and 7d,sim/7max,sim ~ 0.0098 
for q — 10^" model, respectively, i.e., once the kinetic energy of particles become 1% of gt^imax, the radiation drag force is 
comparable to the Lorentz force. As a result, the radiation drag force for the super particle would be exaggerated. To perform 
radiation drag force properly, we introduce the reduction factor rj, namely 7^ in (|22[1 replaced by 777^ and take 77 — 0.05 for 
q = 10^^ model and 77 = 0.089 for q = 10^® model respectively although we simulate with 77 = 1 in previous work. 



2.7 A Test Run: Reproduction of quiet solution 

The pair creation in the gap has a significant role in the structure of the magnetosphere. If p air creation is suppressed, the 
static magnetosphere by KM is reproduced. Following up to the Jackson's gedankenexperiment (|jacksonlll97i ). our simulation 



starts with the condition that the positively charged magnetized rotating conductive sphere is initially put at the origin in a 
vacuum. The charge separated particles are emitted from the stellar surface and are accelerated by the induced electric field 
along the magnetic field lines. But the energy of the particles is lost due to the radiation drag force with time. The particles 
are located at the bottom of the potential along the magnetic field line, i.e. the force-free surfaces where the magnetic field 
aligned electric field vanishes. In vacuum, the force-free surfaces are polar domes and the equatorial plane. The inside of the 
polar force-free surfaces are filled with the negative charges emitted from the polar caps of the stellar surface, and they extend 
upward and form north and south force- free domes. Similarly, the positive particles accumulate above and below the equatorial 
plane and form a force-free torus. Eventually the field aligned electric field on the stellar surface is shielded, and the static 
negative charged domes above the poles and the static equatorial positive charged disc are formed. Fig. [T] shows the static 
particle distribution with the equipotentials (j>- For comparison, they are superposed on the co-rotational equipotentials which 
is given by 4>co = 0,fnp/c, where ip = fi sin^ 6/r is the dipole magnetic flux function. The equipotentials follow the co-rotational 
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Figure 1. The particle distribution and the equipotentials for a — 1 and Qunit — 

10~^ model. The orange and cyan points are positive 
and negative particles on the meridional plane (ro, z), respectively. The light cylinder is locates at ro = 5 with the green line. The broken 
lines are the corotational equipotentials. The solid lines are the equipotentials. The contour level is taken with equal interval at 0.5(^efl- 



equipotentials in the negative charged clouds and the inner part of the positive charged torus close to the stellar surface. In 
the cusp of the torus, the two kinds of equipotentials are deviated but they are parallell. Roughly, the E\\ is screened out in 
the both signed charged clouds. Meanwhile, the electric field in the vacuum region faces the direction such that particles draw 
back to the cloud with the same sign, and therefore the cloud-vacuum boundaries are stable. 
The force-free plasma in the cloud has toroidal velocity hy E x B drift, which is given by 

vt = ujft + czu , (25) 
dip 

where zu is the radius in the cylindrical coordinates and (jinco = (p — <j>co- The coincidence of and 4)co in the cloud means 
d0nco/d-(/) — 0, i.e., the co-rotational motion of the plasma. The value of the potentials fall off faster than the co-rotational 
one in the outer part of the disc, and therefore the term cwd(finco/dtp gives super-rotation in the (|25|l . In this structure, the 
toroidal velocity of the plasma in the outer part of the disc was less than ligiit speed because of the size of the disc is much 
less than the light cylinder. 

As mentioned in subsection 12.51 note that our simulation has artificial thin vacuum layer between the inner boundary 
and the stellar surface. The region causes deviation of the toroidal velocity from corotation in the polar dome. Fig. [2] shows 
d4>/dTp for the 10~^ model and the 10~® model in unit of Q/c on the inner boundary. This curve has a symmetry on both 
sides of the equator. Most of the cloud co-rotates except for near the poles and the equator. The maximum deviation from 
the co-rotational velocity is 40% for the 10~^ model and 30% for the 10~® model, respectively. Although the error is large on 
the equatorial plane, an affect near the equator is trivial because the co-rotational velocity is originally small at this place. 
But the deviation by the numerical error near the pole is significant if polar electric domes expand in the vicinity of the light 
cylinder. In our simulation, the polar force-free electric dome, which connects to the stellar surface along the magnetic field 
anchored at 3 deg for 10~^ model and 1 deg for 10~^ model, has about 10% sub-rotational velocity from the light speed with 
the distance of light radius. 



2.8 The treatment of pair creation and initial condition 

Fig. [3] shows the distribution for strength of the electric field along the magnetic field line ) on the meridional plane for 
the static solution. A strong intensity region is appeared in the middle latitudes, while the less intensity is appeared in the 
charge clouds. For the static solution, there is a vacuum gap in the middle latitude with a potential drop of ~ 2(/)cff. In the 
vacuum gap, the maximum intensity of is typically -Ey^max ~ 3-Bi at {r^,z) — {1.8R,1AR) on meridional plane, where 
Bi = fi/Rf is the dipole magnetic field intensity at the light radius on the equator. If the charge particles are injected into 
the gap, they will be immediately accelerated to ultra-relativistic regime by the E^^, and radiate the 7-rays by the curvature 



8 T. Wada and S. Shibata 




0.5 



2 4 6 8 10 



10 20 30 40 50 50 70 80 90 

theta 



Figure 2. d<f>/dip nomalized by Q/c versus the colatitude on the inner boundary. The small panel is close up near the polar region from 
Odeg to lOdeg. 
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Figure 3. The Ey intensity map for the pair starved static structure with q = 
intensity is nomalized by Bi. The light cylinder locates at ta = 5 with green line. 



model on the meridional plane, where the color 



radiation process, whose power is given by 



(26) 



37?? ■ 

The typical Lorentz factor of the accelerated particles is estimated by the force balance between the radiation drag force and 
the electric force, i.e., qE^^ = Ec/c, which gives 

7sat = [ ■ (27) 



\2qR 



Then, the typical energy of the emitted 7-rays is given by 



sat 



With typical parameters of the pulsars, we obtain the Lorentz factor of 7aat = 1.1 X 10'^Po//(i?c/iii)'/^(-B||/Bi)''''' and the 
photon energy of — 1.4 x Pq(^ {Rc/ Ri)^^'^ {E^^ / Bi)^^* GeV. In the present work, we simplify the pair-creation process, that 
is, (1) we do not distinguish between the photon-photon pair-creation process and the magnetic pair-creation process, and 
(2) we ignore the effects of the collision angle between the 7-rays and soft-photons (or the magnetic field). In our simulation, 
we determine the pair-creation position using the condition that > 2m(? , implying the pair-creation is occurred at the 
position, where the electric field is stronger than E^v = 6.6 x 1Q~^ Pq2{.Rc / Rif^^^ Bi. e-y > 2mc^ is not applicable in reality 
because of the effect of collision angle of the two photons. On the other hand, we shall see below that actual setting of the 
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name q Ea Tpr Ar A6 



AO 


10-^ 


10-3 


2.0 


0.1 


2° 


Al 


10-^ 


10-3 


1.0 


0.1 


2° 


A2 


10-^ 


10-3 


0.2 


0.1 


2° 


A3 


10-^ 


10-3 


0.1 


0.1 


2° 


A4 


10-^ 


5 X 10-* 


2.0 


0.1 


2° 


A5 


10"^ 


2.5 X 10-" 


2.0 


0.1 


2° 


BO 


10-6 


10-3 


0.2 


0.1 


1° 



Table 1. The numerical parameter of simulation. 



value of Ed is much larger than this value in the numerical simulation since the accuraly of the elelctric field is much larger 
than Ed- 

For our simulation, the pair creation is performed in the following way. We introduce the critical electric field intensity to 
generate pairs in the simulation as parameters, that is Ed/ Bi = 0.25, 0.125, 0.0625 although they are larger than the realistic 
value. At first, the equally-spaced grid points {fi,Oj) are prepared with spherical coordinates in the meridional plane, where 
Af = 0.1, Ae = 2° for q = 10-^ model and Af = 0.1, Ad = 1° for q = IQ-*^ model, respectively. If the field-aligned electric 
field at a grid points E^^{i,j) is larger than Ed, then we put on n± = njvi pairs around the grid, where um is multiplicity of 
the pair creation. The generation of pairs is repeated in every -fpr, where we take fpr — 2.0, 1.0,0.2, which are less than the 
light crossing time of the light radius. We consider um as a constant parameter and the lowest cases um ~ 1 are carried out 
to guarantee the pure dipole magnetic field assumption. We set up the static solution as an initial condition of the simulation 
with pair creation effect, e.g., 3000 pairs are generated in the gaps with q = 10-^ model with um = 1 in initial. Table [1] 
provides the numerical set of parameters of simulation. 

Note that the perturbation of a simulation particle makes error electric field. The maximum error of electric field caused 
by a simulation particle is roug hly given by Eerr = g/A^ which is 0.25Bi for 10-^ model and 0.025Bi for 10'^ model in the 
present parameter setting, where A is typical interval of the grid to estimate and we choose A — O.IR. If Ed is taken to 
be much smaller than -Eerr, the pair is generated excessively and therefore the calculation is immediately broken by running 
out of the limit of the number of particles in our simulation. 



3 RESULT 

3.1 Quiet Solution: without pair creation 

As has been shown, if pair creation is suppressed, the electrosphere is composed of electronic domes above the pole, a positronic 
or ionic equatorial disc and vacuum gaps in the middle latitudes. It takes about one rotation period until the quiet state 
is achieved with a — 1. We also simulate for smaller system charges of a = 0.75,0.5,0.4,0.3 with q = 10~^ model. The 
smaller the net charge, the larger the extent of the polar dome (see left panel of Fig. |4|. This tendency has been stated by 



Krause-Polstorff fc Michell (|l985bl ) + 10 and -f4 models which correspond to our a = 1 and 0.4 models, respectively. 

The same structures are identified as well as KM with a ^ 0.4 in our simulation because the polar dome falls in the 
light cylinder and therefore it satisfies the condition of rigid constraint for the particle by the dipole magnetic field. But the 
model with a — 0.3, the static structure is changed. The global flow pattern of the particles for the model with a — 0.3 is 
shown in the right panel of Fig. |4l When we choose a model with a < 0.4, the polar dome extends beyond the light cylinder 
and the azimuthal velocity of the particle in the dome increases with the axial distance although the azimuthal velocity is 
slightly small from the corotational velocity. Thus, the Lorentz factor increases close to 7d near the light cylinder and then 
the radiation drag force of the particle makes f^^^ x B drift just out side of the light cylinder. Thus, the particles emitted 
along the polar magnetic field line migrate into the inner magnetic surface if they pass through the light cylinder, and they 
returns star just higher colatitude of their departure point. As a result, there are outfiows of the negative particle from the 
polar region, which are bounded the magnetic surface footed on the stellar surface with 4.5 degree. Note that, for our case, 
the whole charge cloud is still in magnetic field dominant region {B > E). When the outflows go beyond the light cylinder, 
they move across the inner magnetic surface by the f^^^ x B drift and they become inflow by the quadrepole electric field 
in the vacuum region at the middle latitudes. Thus, the closed poloidal current loops above the pole are formed in several 
light radius. This is more fav orable rather t han the faraway loop of the electron made by attraction of the monopole electric 



field of the star suggested by Jackson ( 19761 ). Our result confirmed the same flow pattern of electrons given by Rylov (1977) 



Although a part of the edge of the dome has fast azimuthal velocity, the current density of the dome in the vicinity of the 
light cylinder is typically 10-''pGJ,poicC, and therefore the modification of the dipole magnetic field should be ignored. 

Fig. [5] shows (j) along two radial directions with 6 = 30° and = 90° (equatorial plane) , respectively. The former radial 
line passes through the polar dome, while the latter does through the equatorial disc. It is notable that the derivative d^/dt/" 
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Figure 4. Left panel: The distributions of the particles on the meridional plane with simulated Qsys's, Colored dots indicate a = 
0.75, 0.5, 0.4, 0.3 with blue, magenta, orange and yellow respectively. The light radius is located at ro = 5 with green line. Right panel: 
the current density pattern on meridional plane for a = 0.3 model. The arrows indicate the direction of the current density of negative 
charged particle(blue) and positive charged particle(red). The magenta curve indicates magnetic field line footed on the surface with 
colatitude 4.5 degree. The boundary surface with E = B is drawn by solid curve. 




indicates the angular velocity of the clouds. Along the line with = 30°, the electric potential (j>{tp) follows the co-rotational 
<l)co in between A and B in the figure, which correspond to the stellar surface and the surface of the cloud (cloud-vacuume 
boundary). This indicates that the cloud co- rotates with the star in the cloud. Along the line with 6 — 90° (equator), the 
inner part of the disc between D and E follows the co-rotational (jico, and therefore they show co-rotation and E\\ = 0. The 
point F corresponds to the top of the equatorial disc. Because of d(j)/d^ > Q/c in between E and F, the positive charge at the 
edge of the disc are in super-corotation. This part is connected to the vacuum gaps in the middle latitudes along the magnetic 
field line. The electric potential in the outer part of the disc decreases faster than the co-rotational potential. Thus, the super 
corotation of the disc top is obtained. Outside of the clouds (between B and C) is in vacuum, and deviates from 4)co , i-e., 
_E|I exists. These features can actually be seen with the velocity of the particle in Fig. [6l i.e., one can see the co-rotational 
motion in the polar domes and in the inner part of the disc, the super-rotation in the outer part of the disc. The azimuthal 
velocities are not close to the light speed at the top of the equatorial disc, and therefore the radial diffusion mechanism by 
rotational inertia and radiation drag are negligible in such cases. Outside of the disc (between F and G) is in vacuum. 



3.2 Active Solution: with pair creation 

The steady solutions are obtained in about 4 rotation periods. For comprehending the trajectories of the particle, we have to 
know the magnitude relation between electric field and magnetic field and the effect of the drift motion caused by radiation 
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Figure 6. The particles color-coded with toroidal velocity on the meridional plane. The meaning of the colors is shown in the color 
bar just below the equatorial plane, where each color indicates the co-rotation speed at each axial distance. The solid curves are dipole 
magnetic field lines. 



drag force. When the plasma is in the magnetic field dominant region (MDR), it tends to move along magnetic surface. 
However, once it is in the electric field dominant region (EDR), particles are accelerated toward the direction of tire electric 
field, ff -Ex ---^ _B in MDR, the velocity of the E x B drift of the plasma is close to light speed. Then the radiation drag force 
for the particle is comparable to the Lorentz force, that is the Lorentz factor of the plasma becomes 7d, which is given by 
(|28|l . Such plasma drifts perpendicular to the magnetic surface by the f,.^ x B drift. It is notable that the f^^^ x B drift is 
opposite direction depending the sign of the charge, i.e., the positive charged particle drifts outward and the negative charged 
particle drifts inward. 

In the left panel of Fig. [T] the typical trajectories of particle on the meridional plane are drawn with red and blue curves 
and EDR and MDR are bounded with thick curves. The EDR appears around the equatorial plane in a wedge-like shape 
beyond uj ~ 4.5 — 0.9Ri with the opening angle of about 50° from the equ atorial plane. Although the EDR appears due 
to the monopole electric field in the quiet model ( Jackson 19761 : Rvlov 1977 ). our present model has almost no net charge. 



The EDR seemed to have been formed by the global structure of the charge clould, and in particular it would be due to the 
growing eq uator ial positive charged disc. The similar shape of EDR is discussed in the force-free model around Y-point by 
Uzdensky The angle between the EDR and the equatorial plane in our result is just wider than Their EDR. It is 

interesting that the structure of EDR is very similar although our result precludes the force-free approximation and assumes 
the magnetic field to be dipole elsewhere. In the right panel of Fig. [71 the luminous color-coded particle have the Lorentz 
factor being comparable with 7d, i.e., the radiation drag force for the particle is comparable to the Lorentz force qBi. They 
are mainly in the vicinity of the light cylinder. The positive charged particle near the equatorial plane just inside the light 
cylinder has the Lorentz factor with 7d. At the same time, the negative charged particle in the vicinity of the light cylinder 
with the height above lOR. The dusk color-coded dots in the right panel of Fig. [7| indicates the non-accelerated particles, 
i.e., they are co-rotating equatorial disc and conically-shaped polar domes. The background color of the left panel of Fig. [7] 
indicate the intensity of E\\, which are color-coded red {E\i > 0), blue (Sy < 0) and white (E^ = 0), respectively. Both sides 
of the outer gap, the E\\ is shielded and there is the co-rotating charged clouds. 

In the left panel of Fig. [71 the flow of the positive particle generated in the outer gap goes into EDR, it easily goes out 
beyond the light cylinder (ao — >■ fli). The plasma at the outer edge of the disc at bo has a fast azimuthal velocity, and therefore 
they are slipped out by f^^^ x B drift (foo — >■ fei). Meanwhile the flow of the negative particles generated in the outer gaps 
moves back to the star (co — ^ ci) with colatitude 24° < < 34° on the stellar surface, and re-emitted from the polar region. 
The polar flow is separated by the magnetic surface footed on the stellar surface with colatitude 10° , where the equipotential 
surface of (f) — connects up to the star. For the lower colatitude region, it is outflow drawn by do — di . For the other region, 
it is circulation on the meridional plane, which is starting from the polar annulus with the colatitude 10° < 9 < 12°, and 
returning to the annulus with colatitude 18° < 9 < 24° (eo — s> ei — 62 — >■ 63 — > 64). On the way through eo — > ei, the flow 
having fast azimuthal velocity migrates to the inner magnetic surface because of f^^^ x B drift. Once it moves in EDR, it is 
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Figure 7. Left panel: typical trajectories of the particle and _E|| intensity map, the curves are color-coded with blue (for the negative 
charged particle), red (for the positive charged particle) and the color map indicates electric field intensity normalized by Bi. The thick 
solid curve indicates isosurface with E = B, and the solid curve with green means equipotential surface with cj> = Q. Right panel: Lorcntz 
factor of the particles, which are color-coded with the Lorentz factor normalized by 7(j. For both panels, the light cylinder is located at 
ro = 5 and the line in the middle latitude is null surface. ) 



accelerated and passes over the equatorial plane. Because of the plane-symmetric fashion about the equatorial plane, the flow 
having the same trajectory in the opposite hemisphere returns from e2, then it moves in MDR (e2 — > 63). Thus, it returns 
along a magnetic field line just outside the pair creating region in the middle latitude (63 ^> 64). 

The overhead view of the steady solution is shown in Fig. [S] The intensity ratio of the dipole magnetic field and the 
magnetic field made of the magnetospheric current in the steady state is drawn with gray scale map. The current density 
normalized by Q,B/{2ti) on the meridian plane is drawn with color-coded arrows and the electrostatic equipotentials normalized 
by the 0eff are drawn with the curves. The green curves indicates that the value of the equipotential is zero. The value of 
the potential is negative in the lower colatitude region and is positive in the higher colatitude region. The contour level is 
incremented by —3, —2, —1,0 with common logarithm from the zero-equipotential surface for both signs of potentials. It shows 
the assumption of dipole magnetic field is valid in the light cylinder and especially round about the equatorial plane, although 
the magnetic field made by the magnetospheric current is comparable to the dipole field in the middle latitude about beyond 
3-Ri, where magnetic fiux is modified to be opened. The circulation pattern of the fiow in the present result would turn down 
if the modification of the magnetic field is concerned. The magnetic field made of the magnetospheric current is trivial inside 
the light cylinder, namely the structure of the outer gaps and pair creation rate should not be affected in present model if the 
modification of the dipole magnetic field is considered. In contrast, outside of several light radius, the magnetic field made of 
the magnetospheric current is comparable to the dipole magnetic field. However, there are almost in EDR, and therefore the 
particle flow is controled by the electric field. Concerning the global structure of the current and equipotential surfaces, the 
negative charged flowed out in a lower latitude and the positive charged flowed out in a higher latitude. However, the system 
of the charge of our result is almost neutral, that is the monopole electric field would not prevent reversed sign outflow at a 
great distance. Eventually both outflows have enough kinetic energy larger than the potential energy at the outer boundary, 
and therefore they can reach inflnity beyond the outer boundary. Inside of the sphere about with radius 3i?i, there are poloidal 
current loops caused by f^^^ x B drift from equatorial plane to the pole, which mainly consist of negative charged flow. In 
previous paper, the parameter of the radiation drag force, 77 is taken to be unit so that the circulation of the both signed 
charge caused by the f-^^^ x B drift is incident. In present result, there are few returning positive charges although the poloidal 
current structure is similar to that in previous work. 
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Figure 8. The arrows indicate the current density made of positive charge (orange) and negative charge (cyan), which are normaUzed 
by f!i3/(27r). The intensity ratio of dipole magnetic field and the magnetic field made of the magnetospheric current is drawn by gray 
scale color-coded map, that is The curves are electrostatic equipotentials normalized by (/>max with common logarithm, where 

(/) = is colored with green, the solid curves indicates <^ > and the dashed curve indicates <^ < 0. 



4 DISCUSSION 



4.1 Comparison with the quiet solution 

Smith. Michel, fc Thackeil (|200ll ) concluded that their result for the axisymmetric magnetosphere is quiet because the pair 



plasma which moves along the dipole magnetic field lines screens out the E\\ in the outer gaps a nd therefore t he pair creation 
is stopped. They implied that activity of pulsar is essentially caused by obliqueness (see also iMichell 120041 ). Although the 
similar screening of the outer gaps is shown in our simulation, but intensity of the electric field is maintained to generate pair 
plasma. We note that the our act ive solution arise from taking smaller threshold intensity of electric field for pair creation than 
Smith. Michel, fc Thacken (|200ll ) and solving the equation of motion for each particle without the assumption of restraint of 



plasma along the dipole magnetic field line. As a result, we showed that if the radiation drag of the particle and the supply 
of the pair plasma in the gaps are performed, the quiet system should be broken and migrates over to the active system with 
gaps and both signed outfiows of the plasmas nevertheless the rate of the pair creation is taken to be small as the modification 
of the dipole magnetic field is omitted. Some fraction of the positive charge generated in the outer gaps accumulates the 
equatorial plane, and therefore the super rotating disc is growing and the azimuthal velocity increases, and then the Lorentz 
factor increases as the light cylinder reaches in some fraction of 7d. The particles at the cusp of the disc start to leak out due 
to the X S drift. If the potential drop in the outer gap is / in fraction of the effective voltage, the trans-field drift would 
take place beyond (1 — /)-/?!, i.e., it is likely that radiation drag causes trans-field drift motion within the light cylinder. As 
has been shown, this effect creates an outfiow of positive charged particles even if the dipole magnetic field is closed, and it 
is favorable to keep the outer gaps. 



As diocot r on in stability of differential rotation disc is pointed by ISpitkovskvl (120041 ): iPetri. Hevvaerts. fc Bonazzola 



( 2002al ) : iPetril (|2007l ). the leaking of the edge of disc is realized after several rotation periods for our simulation. Although it 
prompts decrease of the charge from the system even if all particles are trapped in closed magnetic surface or pair creation is 
suppressed, but the speed of the growth of the disc decreases when the disc grows in vicinity of the light cylinders and the 
decrease of charge is trivial compared with total charge of the disc. We simulated to confirm how the diocotron instability 
affects the global structure of the disc over 20 rotation periods for the pair starved electrosphere, which is barely acceptable to 
carry out for our present environment. At least we realized that the diffusion by the diocotron instability was not so significant 
in such a time scale, although the diffusion by f x B changes the quiet electrosphere in several rotation periods. Note that 
the growth of the disc is made by the pair created positive charges in the outer gaps and the main component of leaking of 
the disc is made by f x S drift of the particle in our simulation and even the loss of the positive charge from the edge of 
the disc is much lower than the outfiow from the outer gaps to EDR (See, left panel of Fig. [71 ao a\). 

If copious plasma is supplied to the equatorial disc from the outer gaps, another important issue remai ns the region, th e 
so-called Y-point, is expected to have the electric field larger than the magnetic field in the force free theory ( Uzdenskvll2003l ) . 
In the next paper, we treat higher rates of pair creation, for which modification of the magnetic field becomes significant in the 
vicinity of the light cylinder. Then, some poloidal magnetic flux is opened and toroidal magnetic field would be anti-parallel 
on both sides of the equatorial plane in which the dissipation of the magnetic field is expected to accelerate the particles by 
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dissipation process of the magnetic field. Thus, trans-field leakage of the particle around the Y-point is a very interesting 
issue, but we postpone further discussions until the simulations with higher rates of pair creation are performed. 



4.2 Formation process of steady solution 

Although the steady solution is stated in section 13.21 here we discuss by the stage of the active solution. Once pair creation 
is on set, the pair plasmas are generated at first in middle latitudes in which _E|| is stronger than Ec^ (see Fig. [3]). The pairs 
created in the outer gaps are immediately separated in oppsite directions along a magnetic field line by E\\ . Then, most of the 
positive charges move into EDR to form an outflow. At the same time, the negative charges move inward. They returns to the 
star and are re-emitted from the polar regions so that the height of the dome increases. The outflow of the positive particles 
causes significant change to the global structure: the system charge, which is initially assumed to be positive, is reduced. As a 
result, the potential of the polar region becomes negative, and a part with lower latitude of the polar domes becomes outflow. 
The other half of the polar dome grows and crosses over the light cylinder. The Lorentz factor of the plasma increases to bring 
about trans-field motion by radiation drag force. It is notable that the negative charges obtain kinetic energy if they drift 
toward inner magnetic surfaces. Some part of the particles become outflow if the kinetic energy is larger than the potential 
energy at the outer boundary, and therefore the separatrix of the negative charged flow is made on the equator with 3_Ri 
in present result although it is set on the equator about with 8-Ri in previous work. The difference might be made by the 
resolution of the charge and the detailed mechanism of the formation of the separatrix with modification of magnetic field 
needs additional simulation in future work. Thus, steady state is found in equal losses of both particle species and the lost 
particles are supplied by the outer gaps. 

In our simulation, we demonstrate the generation of the particles intermittently in some periods fpr = 0.1,0.2, 1.0,2.0 
for the pair creation and a fixed period fpo = 0.1 for the popping from the stellar surface. All steady results had the same 
geometry of the charge clouds, which were equatorial disc of positive charge, polar domes of negative charge and both signed 
outfiows. Thus, the choice of -fpr and fpo in present simulation does not affect the results. We checked that if the lower system 
charge solution without pair creation is chosen as a initial condition, then the same result was obtained qualitatively. 



4.3 Structure of Outer Gap 

We here discuss the structure of the outer gaps for our results. For the discussion, the structure of the electrostatic potential 
and the intensity of E\\ in the pair creating region should be realized for our active solutions, where we obtained the solution 
with the several rates of the pair creation in simulation. Note that in our simulation, the Ecr is taken to be large artificially 
compared with a parameter for emission of gamma-ray to reality, in other words, pair creation rate is underestimated. For 
the reason, the size of the pair creating region is expanded and therefore the geometry of the outer gaps is not revealed 
quantitatively from our simulation. Although the structure of the outer gaps which are demonstrated the lower rate of the 
pair creation have contained some fraction of artifact, the result implies how the transition of the magnetosphere from the 
active state to the quiet state is realized with time for the pulsar, that is the evolution of the magnetosphere with the 
abundance of the pair plasma, which decreases with age of the pulsar. 

Another important point, we can suggest the possibility of coexistence for polar cap and outer gap by comparing with 
the deviation from the force-free magnetosphere. If a force-free condition are satisfied elsewhere, we give the co-rotational 
potential under dipolc magnetic field assumption, 

0co = ^^^. (29) 
cr 

Then, the non-corotational potential defined (fi-aco = (t> — 4>co- Fig. [9] shows the non-corotational potential with color map 
and pair creating region with contour for Bq model on the meridional plane. The pair creating region become dramatically 
small compared with initial state in our simulation but the intensity of the E\\ in steady state remains larger than Ecr in the 
countour. Eventually, the pair creating regions remained in the middle latitudes and just above the poles. We don't distinguish 
between the process of the pair creations in present simulation, which are photon-photon collision and magnetic pair creation, 
and therefore the rate of the pair creation is only proportional to the product of the intensity of the and the volume of 
the gap. Thus, the supply of the plasma generated in the polar pair creating region was trivial compared with the one from 
the outer gaps, th at is, our simulation underestim ated the effect of shielding by pair plasma with magnetic pair creation in 
the polar cap. As iTakata. Wang, fc ChengI (|2Q10l ') suggested that a new outer gap closure mechanism by the magnetic pair 



creation near the stellar surface is significant to realize the observed features of the gamma-ray pulsars, the effect of the 
different type of the pair creation should be investigated in future work. 

In Fig. [9l the equatorial disc and polar domes has domain in which cj>-nco have deviation less than 0.01(^cff. They are not 
accelerating region of the particle; dead zone. For our result, the existence of the dead zone provides a inspiration. Usually 
conventional polar cap and outer gap are defined on the last closed field line for GJ model and the poloidal current in the 
gaps directs oppositely, and therefore the gaps can not coexist. In present result, the dead zone in equatorial disc is reduced 
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Figure 9. Distribution of non-corotational potential on meridional plane, and contour indicates isoline with _E|| 
is normalized by 2% of the open field line voltage 4>cS- The line in the middle latitude is null surface. 
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in the light cylinder and the effective last closed line is migrated to higher colatitude and the one in polar dome cut down the 
radius of the cone in which particle outflows. Note that the dead zones divide the last closed field line and would be able to 
coexist polar cap and outer gap on different magnetic surfaces. The poloidal current travels into the polar region of the stellar 
surface of the lower colatitude side of the polar dead zone, migrates to the magnetic surface with higher colatitude, and passes 
through the outer gaps. The current makes poloidal circulation beyond the light cylinder. Because we do not demonstrate 
detailed processes of pair creation and do not simulate with enough accuracy to resolve the size of the polar caps in present 
results and therefore it is not clear whether current from the polar cap is necessary condition or not. This interesting result 
should be discussed in future work. 

For confinement for the robustness of our gap-wind solution, we carried out the models of some rates of the pair creation. 
We checked that E\\ in the gap saturates with a constant value for all steady states. Fig. 1 101 shows the saturated intensities of 
on the null surface in all simulations. All results had regions whose E\\ was larger than Ecr on the null surface. The width 
of the base of the bell-curve on the line with E\\/Eci = 1 indicates typical size of the pair creating region on null surface. 
There is a tendency for the higher pair creation rate model {Ao — >■ Ai — >■ A2 — >■ A3), which is parameterized by the frequency 
of the pair creation tpy, to decrease the size of the pair creating region and the intensity of E\\. The saturated intensities 
are maintained just above the E^- When E^t is taken to be much smaller value (Ao — >■ ^45), the pair creating region 

decreased in the same way and the saturated electric field intensity was held just higher than E^- In the lower E^t model, 
the particles tended to screen out the gaps near the star although the outer edge of the pair creating region does not change. 
This indicates that the inner boundary of the outer gap, which a has beak-like shape on the meridional plane, would not be 
connected to the star if we could carry out a much lower Ea model. Considering from figure [71 the inner edge of the outer gap 
would be refined just outside the outer edge of the dead zone which is roughly defined as a magnetic field line footed on the 
stellar surface with 32° in our simulation. The inner boundary is slightly ins i de co mpared with that of the conventional outer 
gap model on the null surface. This agrees with Takata. Shibata. fc Hirotanil (2004) which is calculated the electrodynamics of 
an outer gap on the meridional plane with assumed external current, then as the current density increases, the inner boundary 
of the outer gap shifts toward the stellar surface. As previously explained, the main path of equatorial outfiow exists outside 
of the dead zone, outer gap structure in our simulation can be maintained steadily. 

We emphasize that the thickness of the pair creating region transversing the dipole magnetic surface would be related 
to the size of the region in which leakage of the disc particles takes place (~ fRi) because the field aligned potential drop in 
the outer gap contro ls the super-rota t ion o f the disc. Phenomenologically-postulated, the thickness of the outer gaps for old 
pulsars are shown bv lZhang fc Chenelj 19971 ) in which they pointed out the effect of the finite mean free path of the gamma-ray 



for pair creation is important for defining the thickness of outer gap. In addition to this, the direction of the gamma-ray would 
be important to define the geometry of the outer gap because the particle has fast azimuthal motion in the gaps near the 
light cylinder, and therefore the gamma-ray is emitted to the azimuthal direction and the path of the gamma-ray should be 
considered with three dimensional geometry. We should consider the mean free path of gamma-ray for pair creation with a 
three-dimensional model in future work. 
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Figure 10. The curves are normalized by Ecr on the null surface. The horizontal axis indicates distance normalized by the stellar 
radius along the null surface. The tags indicate type of simulation referred in Table [T] 



4.4 Polar Gap 

The polar cap accelerator has been an outstanding issue in the context of how one can understand radio emissions. We find a 
potential drop just above the stellar surface in our simulation in a similar fashion to the WS. The mechanism is very simple: in 
a steady state the negative charges created in the outer gap flow back to the star, and therefore the equal amount of negative 
charges should be emitted for steadiness. Then, he potential drop in the polar region of the stellar surface controls the polar 
flow of negative particles, which is represented by 

In our all results of simulations, <j>poic becomes slightly negative, and therefore the negative particles are pushed and easily 
escape to infinity. Implied by (|30p . this potential is kept by balance of the numbers of negative and positive particles in the 
magnetosphere (the second term) and Qsys (the third term). Thus, we see that in the steady state Qsys is determined so that 
the losses of negative and positive particles are balanced. This indicates that the polar caps have a close correlation to the 
injection of the electron from the the outer gaps. In the present case, the rate of the emitted particles from the polar region is 
less than would be capable of shielding the electric field above the pole, that is the current density from the pole is less than 
Pgjc, and therefore the unscreened potential increases gradually along polar magnetic field lines. In other words, this potential 
drop above the pole is not a conventional polar cap because the electric field is not shielded in finite distance near the stellar 
surface. There is a very interesting issue whether a potential drop near the polar region remains if enough electrons return 
from the outer gaps and are re-emit ted from the pole or genera ted in the polar cap and then the emitted particles should 
form the space charge limited flow ( Arons fc Scharlemann 19791 ). but our simulation particle has large inertia length with 



3% of stellar radius, and therefore the detailed structure might not be demonstrated. Note that our simulation contains the 
artiflcial polar gap arising from the intermittent emission of particle on the stellar surface, which is in the range of fpoC = O.li? 
at a maximum. Unfortunately, it is difficult to take a much smaller value of fpo in the present parameter setting. Thus, it is 
not clear that polar electric field is shielded on the pole if pair creat ion process is consid ered in detail and whether the polar 



electric fiow needs pair creation in polar cap or not. As shown by ( Hirotani et al. 20031 ). the current from polar cap have a 
significant effect on outer gap electrodynamics, to discuss the possibility of existence of the polar caps linked with the outer 
gaps, the combination of localized simulation and global simulation would be needed. 



5 CONCLUSION 

For an axisymmetric pulsar magnetosphere, pair creation in the outer gaps results in expansion of rotating electrosphere 
and the trans-field drift motion due to radiation drag force near and beyond the light cylinder. Eventually, a steady state is 
achieved with the outflow of both particle species and global current loop on the meridional plane. We conflrmed the pair 
starved static electrosphere shifts to the steady structure by the pair plasma in the gaps. The pair creation plays important 
role of metamorphosing from the static electrosphere to the steady structure nevertheless the pair creation rate is artiflcially 
suppressed in our present model. For more detailed discussion of the structure of the magnetosphere, the consideration of 
the mean free path of the pair creation, which is important for old pulsars having thick outer gaps, is needed to discuss the 
structure of the gap quantitatively. Additionally, a higher pair creation rate case should be simulated, then modiflcation of 
the magnetic fleld from the dipolar changes the pattern of the outflow of the plasma and the structure of the gaps. Thus, our 
model can be developed to consider pulsar magnetosphere including much more pair plasma in the next paper. 
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